#!/usr/bin/perl 
##############################################################
# Program to generate a time series
#
# - parameters for time series are read from a parameter file
#   described below. The name of the parameter file is entered
#   as the first command line option.
#
# - usage:
#      ./dataGen.pl <param file>
#
##############################################################
# Input parameter file: consists of an asci text file with a
#   list of parameters followed optionally by a comment
#   starting with a '#' character.
#   - The format of the file is as follows:
#
#    <length>          # length of time series
#    <start> <stop>    # time range
#    <E|R>             # 'E'ven or 'R'andom sampling
#    <% random drops>  # *if* random sampling, % to drop
#    <# components>    # number of components in series
#    <component n   = (freq amp phase) >
#    <component n+1 = (freq amp phase) >
#
#
##############################################################
# ben.mort@gmail.com                     Last update:  4/09/05
##############################################################
use warnings;
use strict;
use dataGenSubs;    # functions for dataGen

my $pi = 3.141592654;
my @amplitudes;     # amplitude array for time series
my @args;           # sine args used in calculation
my @sines;          # sine wave amp components used in calc.

dataGenSubs::welcomeMessage();

# Check if there are the correct number of command line args.
dataGenSubs::checkNumCmdArgs($#ARGV);

# Process parameter file
#========================
print "PARM:\n";
print STDOUT "PARM:\tGetting inputs from parameter file\n";
print STDOUT "PARM:\t==================================\nPARM:\n";

# open the paramter file
my $parmFile = $ARGV[0];
if ( -e $parmFile ) { print STDOUT "PARM:\tParameter file\t= '$parmFile'\n"; }
else { die "\nERROR:\tparameter file '$parmFile' does not exist!\n\n"; }
open( PFILE, "$parmFile" );

# get the length of the timeseries
my $length = <PFILE>;
if ( $length =~ /^\#/s ) { $length = <PFILE> }
$length =~ /^(\d+)/s
  || die "ERROR:\tLENGTH paramter not recognised\n\n";
$length = $1;
print "PARM:\tLength\t\t= $length\n";

# get range of the time series
my $range = <PFILE>;
if ( $range =~ /^\#/s ) { $range = <PFILE> }
$range =~ /^(\d+\.\d+|\d+|-\d+\.\d+|\d+)\s+(\d+\.\d+|\d+|-\d+\.\d+|\d+)/s
  || die "ERROR:\tRANGE parameter not recognised!\n\n";
my $tmin = $1;
my $tmax = $2;
print "PARM:\tTime min.\t= $tmin\n";
print "PARM:\tTime max.\t= $tmax\n";
$range = $tmax - $tmin;

# get type of time series (spacing type ect.)
my $type = <PFILE>;
if ( $type =~ /^\#/s ) { $type = <PFILE> }
$type =~ /^(RN|E|D|N|R|SF)/s
  || die "ERROR:\tTYPE not recognised!\n\n";
$type = $1;
print "PARM:\tType\t\t= $type\n";

# random drop type parameters (if the type = R)
my $rdrops;
my $nrands;
if ( $type =~ "R" ) {
	$rdrops = <PFILE>;
	if ( $rdrops =~ /^\#/s ) { $rdrops = <PFILE> }
	$rdrops =~ /^(\d+)/s
	  || die "ERROR:\t\% DROPS not recognised!\n\n";
	$rdrops = $1;
	print "PARM:\t% drop  = $rdrops\n";
	$nrands = int $length * $rdrops / 100;
	print STDOUT "PARM:\tno. of drops   = $nrands\n";
}

# Spacing file parameters (if the type = SF)
my @times;
my $nSpaces = 0;
if ( $type =~ "SF" ) {
	my $spaceFile = <PFILE>;
	while ( $spaceFile =~ /^\#/s || $spaceFile =~ /^\s+\#/s ) {
		$spaceFile = <PFILE>;
	}
	$spaceFile =~ /^(\w+-\w+.\w+.\w+|\w+.\w+.\w+|\w+.\w+)/s
	  || die "oh no! an error!";
	$spaceFile = $1;
	if ( -e $spaceFile ) {
		print STDOUT "PARM:\tSpacing file = '$spaceFile'\n";
	}
	else {
		die "\nERROR:\tSpacing file '$spaceFile' does not exist!\n\n";
	}
	open( SFILE, "$spaceFile" );
	my @rawDATA = <SFILE>;
	close(SFILE);
	my $i = 0;
	foreach my $line (@rawDATA) {
		print $line;
		if ( length($line) > 1 ) {    ## slightly dodgy catch for empty lines
			$line =~ /^(\d+\.\d+|\d+|-\d+\.\d+|-\d+)/s
			  || die "ERROR: cannot read data line!\n\n";
			$times[$i] = $1;
			$i++;
		}
	}
	$nSpaces = $#times + 1;
	$length  = $nSpaces;              ## OVER-RIDE LENGTH !!
	print STDOUT "No. of spacing points = $nSpaces\n";
}    # if ( $type = "SF" ) -- spacings file #

# Duty cyle parameters (if type = "D")
my $dlength;
my $dcycles;
my $dspacing;
if ( $type =~ "D" ) {
	$dlength = <PFILE>;
	if ( $dlength =~ /^\#/s ) { $dlength = <PFILE> }
	$dlength =~ /^(\d+\.\d+|\d+)\s+(\d+\.\d+|\d+)/s
	  || die "ERROR:\t\% DUTY parameters not recognised!\n\n";
	$dlength = $1;
	$dcycles = $2;
	$dlength = int $length * $dlength;
	print "PARM:\tdlength = $dlength\n";
	print "PARM:\tdcycles = $dcycles\n";
	my $dtest    = int $dlength * $dcycles;
	my $dpercent = ( ( $dtest / $length ) * 100 );
	my $ddropped = $length - $dtest;
	$dspacing = int( $ddropped / ( $dcycles - 1 ) );
	print "PARM:\tddropped = $ddropped\n";
	print "PARM:\tdspacing = $dspacing\n";
	print "PARM:\tpercent sampled = $dpercent\n";

	if ( $dtest == 0 ) {
		die "ERROR:\tno measurements made with this duty settings";
	}
	if ( $dtest >= $length ) {
		die "ERROR:\t measurements made ($dlength * $dcycles) >",
		  " number of possible time samples ($length)";
	}
}    # ( $type = "D" ) -- duty cycle

# Noise ( type = "N" ) noise added.)
if ( $type =~ "N" ) { die "ERROR: Noise is not yet implimented\n"; }

# Random + noise ( type = "RN" ) random drops and noise
if ( $type =~ "RN" ) { die "ERROR: Random + noise is not yet implimented\n"; }

# get the number of compentents in the time series
my $nComps = <PFILE>;
if ( $nComps =~ /^\#/s ) { $nComps = <PFILE> }
$nComps =~ /^(\d+)/s
  || die "ERROR:\t\# COMPONENTS not recognised!\n\n";
$nComps = $1;
print "PARM:\tNo. components\t= $nComps\n";

# get the spectral components to generate
my @freqs;
my @amps;
my @phase;
for ( my $c = 1 ; $c <= $nComps ; $c++ ) {
	my $tmp_line = <PFILE>;
	if ( $tmp_line =~ /^\#/s ) { $tmp_line = <PFILE> }
	$tmp_line =~
	  /^(\d+|\d+\.\d+)\s+(\d+|\d+\.\d+)\s+(\d+|\d+\.\d+|-\d+|-\d+\.\d+)\s+/s
	  ||
	  ##$tmp_line =~ /^(\d+|\d+\.\d+)\s+(\d+|\d+\.\d+)\s+(\d+|\d+\.\d+)\s+/s ||
	  die "ERROR:\tproblem reading component values!\n\n";
	$amps[$c]  = $1;
	$freqs[$c] = $2;
	$phase[$c] = $3;
	print
	  "PARM:\t  ($c) Amp: $amps[$c]\t  Freq: $freqs[$c]\t  Phase: $phase[$c]\n";
}
print "PARM:\n";
close(PFILE);

# Generate sample time series with even spacings
#================================================
# This step is required prior to generation of several non even spacing types.
# - random drops
# - duty cycle
# - noise/ random noise

if (   $type =~ 'E'
	|| $type =~ 'R'
	|| $type =~ 'D'
	|| $type =~ 'N'
	|| $type =~ 'RN' )
{
	print STDOUT "TIME:\n";
	print STDOUT "TIME:\t==> Generating even sampled time series\n";
	my $tinc     = $range / ( $length - 1 );
	my $tmp_min  = $tmin;
	my $tmp_max  = $tmp_min + $tinc * ( $length - 1 );
	my $freq_max = 1 / ( 2 * $tinc );
	print STDOUT "TIME:\tRange  = $range \t(time range)\n";
	print STDOUT "TIME:\tLength = $length \t(items of data)\n";
	print STDOUT "TIME:\tTime series increment\t= $tinc (range/length)\n";
	print STDOUT "TIME:\tTime series min\t= $tmp_min\n";
	print STDOUT "TIME:\tTime series max\t= $tmp_max\n";
	print STDOUT "TIME:\t- Max. possible freq. in the series = $freq_max\n";


	# generate even time series
	for ( my $j = 0 ; $j < $length ; $j++ ) {
		$times[$j]      = $tmin + $tinc * $j;
		$amplitudes[$j] = 0.0;
		for ( my $c = 1 ; $c <= $nComps ; $c++ ) {
		    $args[$c] = 2.0 * $pi * $freqs[$c] * $times[$j] + $phase[$c] * $pi;
		    $sines[$c][$j] = $amps[$c] * cos( $args[$c] );
		    $amplitudes[$j] += $sines[$c][$j];
		}
	}
}

# Generate sample time series at predetermined times
#====================================================
if ( $type =~ 'SF' ) {
	print STDOUT "TIME:\n";
	print STDOUT "TIME:\t==> Generating spacing file seperation time series\n";
	for ( my $j = 0 ; $j < $length ; $j++ ) {
		$amplitudes[$j] = 0.0;
		for ( my $c = 1 ; $c <= $nComps ; $c++ ) {
			$args[$c] = 2 * $pi * $freqs[$c] * $times[$j] + $phase[$c] * $pi;
			$sines[$c][$j] = $amps[$c] * cos( $args[$c] );
			$amplitudes[$j] += $sines[$c][$j];
		}
	}
}

print STDOUT "TIME:\t==> Done generating time series.\nTIME:\n";




# Generate list of random drops
#===============================
# if random remove $rdorps % points
my @rand_number;
my $n_rands;
if ( $type =~ 'R' ) {
    print STDOUT "TIME:\n";
    print STDOUT "TIME:\t==> Generating spacing random drops\n";
    my $seed = time + int( rand 1000 );
    print STDOUT "TIME:\t-> Generating list of random drop points\n";
    srand $seed;    # seed the random number generator
    $n_rands = $nrands - 1;    # length of random number list to generate
                               # -1 corrects for rand_number array starting
                               # at 0
    my $i = 0;
    print STDOUT "TIME:\t";

    while ( $i <= $n_rands ) {
        $rand_number[$i] = int rand $length;
        for ( my $k = 0 ; $k < $i ; $k++ ) {
            if ( $rand_number[$i] == $rand_number[$k] ) {
                $i--;
            }
        }
        $i++;
    }
    print STDOUT "\n";
}
@rand_number = sort dataGenSubs::by_number @rand_number;






# Write out time series
#=======================
my $tseriesFILE = substr $ARGV[0], 0, -6;
dataGenSubs::checkParamFileType();
if ( $tseriesFILE !~ /dat$/s ) { $tseriesFILE .= ".dat"; }

print STDOUT "DISK:\n";
print STDOUT "DISK:\t==> Writing out time series\n";
print STDOUT "DISK:\tTime series output file: '$tseriesFILE'\n";
open( TFILE, ">$tseriesFILE" )
  || die "\nERROR:\toutput file '$tseriesFILE' creation failed!\n\n";

# Even or spacing time series 
if ( $type =~ 'E' || $type =~ 'SF' ) {
	for ( my $j = 0 ; $j < $length ; $j++ ) {
		my $dline =
		  sprintf( "%.4lf      %.4lf\n", $times[$j], $amplitudes[$j] );
		print TFILE $dline;
	}
}

# Random spacing time series
my @rtimes;    # time series array after random drops
my @ramps;     # amp array after random drops
if ( $type =~ 'R' ) {
	my $drop = 0;
	my $val  = 0;
	for ( my $j = 0 ; $j < $length ; $j++ ) {
		if ( $drop <= $n_rands ) {
			if ( $j == $rand_number[$drop] ) {

				#print TFILE "$j >> drop\n";
				$drop++;
			}
			else {
				print TFILE "$times[$j]  $amplitudes[$j]\n";
				$rtimes[$val] = $times[$j];
				$ramps[$val]  = $amplitudes[$j];
				$val++;
			}
		}
		else {
			print TFILE "$times[$j]  $amplitudes[$j]\n";
		}
	}
}

# Duty cycle time series
if ( $type =~ 'D' ) {
	my @darray;
	my @dIncsArray;
	my $iduty  = 1;
	my $ispace = 0;
	my $Ninc   = 0;
	$dIncsArray[0] = 0;

	for ( my $j = 0 ; $j < $length ; $j++ ) {
		my $incTest = ( $iduty * $dlength - 1 ) + ( $ispace * $dspacing );

		if ( $j == $incTest ) {
			if ( $Ninc / 2 == ( int $Ninc / 2 ) ) {    ## even
				$ispace += 1;
				print "space >> $incTest\n";
			}
			else {
				$iduty += 1;
				print "duty >> $incTest\n";
			}
			if ( $Ninc < ( 2 * $dcycles - 1 ) ) {
				$Ninc += 1;
			}
			$dIncsArray[$Ninc] = $incTest;
		}
	}

	for ( my $xx = 0 ; $xx <= $#dIncsArray ; $xx++ ) {
		print "$xx ---  $dIncsArray[$xx]\n";
	}

	my $ii = 1;
	for ( my $j = 0 ; $j < $length ; $j++ ) {
		if ( $j <= $dIncsArray[$ii] && $ii / 2 != int( $ii / 2 ) ) {
			$darray[$j] = 1;
		}
		else {
			$darray[$j] = 0;
		}
		if ( $j == $dIncsArray[$ii] ) {
			$ii += 1;
		}
	}

	for ( my $j = 0 ; $j < $length ; $j++ ) {
		if ( $darray[$j] == 1 ) {
			my $dline =
			  sprintf( "%.4lf  %.4lf\n", $times[$j], $amplitudes[$j] );
			print TFILE $dline;
		}

	}
}    ### duty cycle type

close(TFILE);
print STDOUT "DISK:\t==> Done\nDISK:\n";



# write out header info file
#============================
my $theaderFILE = substr $tseriesFILE, 0, -4;
$theaderFILE .= ".info";
print STDOUT "DISK:\t==> Writing time series header info\n";
print STDOUT "DISK:\tTime series header file: '$theaderFILE'\n";
open( HFILE, ">$theaderFILE" )
  || die "\nERROR:\toutput file '$theaderFILE' creation failed!\n\n";

my $timeStamp = localtime;
print HFILE "**********************************\n";
print HFILE "*    Time series summary info    *\n";
print HFILE "**********************************\n\n";
print HFILE "Datafile: $tseriesFILE (created $timeStamp)\n";
print HFILE "\n";
print HFILE "Parameter table: \n";

if ( $type =~ 'E' ) {    # even spacing
	print HFILE " -even spacing selected\n";
}
if ( $type =~ "R" || $type =~ "RN" ) {    # random spacing
	print HFILE " -random spacing selected\n";
	print HFILE " -% points dropped: $rdrops\n";
	my $dataLeft = $length - $nrands;
	print HFILE " -# points dropped/left: $nrands/$dataLeft\n";
}
print HFILE " -length: $length samples\n";
print HFILE " -range: $tmin -> $tmax\n";
print HFILE " -number of frequency components: $nComps\n";
print HFILE "\nComponents table:\n";
print HFILE "  #  | amp   |  freq  |  phase \n";
for ( my $c = 1 ; $c <= $nComps ; $c++ ) {
	print HFILE "  $c\t$amps[$c]\t$freqs[$c]\t   $phase[$c]\n";

}

# Write out a list of random drop points (if $type == R)
if ( $type =~ 'R' ) {
	my $r = 0;
	print HFILE "\n";
	print HFILE "Random drop points:\n";
	foreach (@rand_number) {
		print HFILE "$_\t";
		if ( $r != 0 && $r % 7 == 0 && $r != ( $#rand_number - 1 ) ) {
			print HFILE " \n";
		}
		$r++;
	}
	print HFILE "\n\n";
}
close(HFILE);
print STDOUT "DISK:\t==> Done\nDISK:\n";


# Finish everything up
#------------------------------------------------------------------
print STDOUT "INFO:\n";
print STDOUT "INFO:\tTime series generation complete. \n\n";



